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ANALYTICAL HAZARD REPRESENTATIONS FOR USE 



IN RELIABILITY, MORTALITY, AND SIMULATION STUDIES 

Part I 

D . P . Gaver * 

M. Acar 



I. INTRODUCTION 

The failure rate function, or hazard function ( hazard 
for short) may be described as the conditional probability of an 
equipment's failing at operating age t, having survived to that 
age. The reliabilities of a variety of electronic and mechanical 
items are conveniently and naturally described in terms of the 
appropriate hazard function, and so is the longevity of human 
beings. The term force of mortality replaces hazard in the 
latter context. 

This paper is devoted to a study of several simple 
analytical representations for hazard functions. These repre- 
sentations are in turn based upon representations of random 
variables having certain required properties, in terms of others 
having familiar distributions — in particular the exponential. 

Similar ideas are due to Tukey (1976) and recently have been examined 
by Parzen (1973). The hazard representations proposed are quite 
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expeditiously used in simulation studies, e.g. of system 
reliability or availability in terms of component lifetimes. 

They may also be used in data analysis studies, in order to 
parsimoniously describe data sets in terms of perturbations of 
convenient and familiar standard distributions. Their use in 
data analysis and simulation is also described in Gaver, Lavenberg, 
and Price (19 76) , and in Gaver and Chu (19 77) . 



II. SYSTEM FAILURE PATTERNS 

It is plausible to think that the time series of failures 
in a system may involve these stages . 

1. Early failures . There may be a relatively large number of 
failure soon after a system is introduced because of design 
defects, production errors, or errors stemming from 
maintenance personnel inexperience. This situation is 
characterized by a hazard function that is initially large, 
but that decreases with time. "Infant mortality" is in 
evidence . 

2. Random Failures . Following the early failure period there 
may be a period during which failures occur at an essentially 
constant rate for a rather prolonged time. During this 
period the hazard function is nearly constant,- so the times 
between failures are close to being exponentially distributed. 
The effect of age or wearout is not yet apparent. 
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3. Wearout Failures . Eventually following the period during 
which a constant hazard is evident there is likely to be 
a period of ever-increasing failure rate caused by wearout 
of system components. 

A graphical representation of a hazard function that 
exhibits the behavior described is given below. Note that it 
has the legendary "bath tub" shape. 



Hazard 



▲ 




time 



Some comments on the above follow: 

A. The term "failure” may refer to an event that is analogous 

to human death, after which the entire system is replaced. 
On the other hand repair or component replacement may occur 
after failure; the system is only repaired, not entirely 
replaced. In the former case, a hazard function of the 
kind depicted in Figure 1 applies to each system event 
("death"); when the system is installed (or is born) that 
hazard operates starting from scratch at t = 0 until 
system failure (human death, for instance) after which a 
similar hazard goes into effect, starting once again from 
zero. In the latter case, in which repair of a component 
occurs, a hazard function like that of Figure 1 applies at 
t = 0, but after the first event ("failure") at t^, a 
repair action is accomplished. The s ame hazard operates 
for t >_ t^ until the next event at > t^, and so 

on. Intermediate situations may be envisioned, in which 

after event n at t the hazard governing system failure 

n+1 starts at t -t,0<t <t. 

n n n n 

B. Although there is reason to assume that hazards somewhat 
like that of Figure 1 occur in general for systems, the 
possibility exists that the system hazard is "bumpy" 
because wearout failures of components or subsystems 
may well occur at intermediate times. 
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C. If the theory is applied to systems with little or no 
wearout propensity , as should be the case when dealing 
with computer software modules, then the hazard function 
may well exhibit the initial falloff of Figure 1 but not 
the rise at later times. In fact, a constant decline as 
bugs are found and removed could be (optimistically) 
anticipated for software. The right-hand side of the bath 
tub vanishes, and the picture is that of a ski slope. 

III. MODELS FOR THE HAZARD FUNCTION 

In this section mathematical models are presented for the 
failure rate or hazard function. Recall that the hazard may be 
defined as follows. 

Definition . Suppose that the time to failure, X, is a random 
variable with distribution function F(x), where F(0) = 0; 
the latter possesses the density function f (x) , f (x) = dF/dx, 
such that for any positive x. 



Then the hazard function , or failure rate at age x, is given by 



x 



F (x) = / f (y ) dy . 
0 



(3.1) 




(3.2) 
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The interpretation of h(x)dx is that it is the conditional 
probability of failure in the interval (x, x + dx) , given 
that there has been no failure up to age x. 

Express the hazard as 



h (x) 



dF/dx 
1 - F (x) 



h (x) dx = i" - - YJ^~) = ~ d{ l°g [1 ~ F (x) ] } , 



it then follows after integration that 

x 

F (x) = 1 - exp[- / h (y ) dy ] . (3.3) 

0 

Thus if the hazard is specified, so is the distribution function, 
and conversely. 

Note that if 



h(x) = X > 0 , 

then 

F (x) = 1 - e Ax , 



0 .1 x 

0 < x , (3.4) 



so a constant hazard function implies the exponential distribution 
of the random variable X, and conversely. 
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Obviously a constant hazard representation does not 
describe the bath tub hazard shape of Figure 1, nor does it 



represent a situation in which hazards decline, possibly because 
design defects or "bugs" are occasionally removed. Here are two 
hazard representations likely to be useful for such purposes. 

(1 ) A Bath tub Model 

Define the random variable Z in terms of X, X being 
exponentially distributed with mean X ^ , as follows: 



where 

a) L(x) is convex in x, L(0) < 1, L (°°) = 1, 

b) R(x) is concave in x, R(0) = 1, R(0) > R(°°) . 

Then the hazard of Z may be made to exhibit a bath tub shape, 
as in Figure 1, by proper choice of the function L and R. 

Example . Suppose 



Z = G (X) = XL ( X) R (X) 



or 



( 3 . 5 ) 



= X<}> (X) , cf> (X) = L (X) R (X) 




a > 0 , 0 < x 



( 3 . 6 ) 




6 > 0 
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Clearly 



z = xL (x ) R(x) 



x 



ax 

1 + ax 



1 ax _ __ x 

I + 6x 1 + ax 1 + 8x 



is a monotonically increasing function of x. Furthermore, 

_3 

choose a large (e.g. a = 10) and 3 small (e.g. 8 = 10 ). 

Then it is intuitively clear that (i) small x-values transform 
into even smaller z-values, e.g. x = 1 corresponds to z = 0.91 
and x = 2 corresponds to z = 1.90, but (ii) this effect 
dwindles as x increases, so x = 10 corresponds to z = 9.8 
and x = 50 to z = 47.5 and the z-values closely resemble 
the x's percentage-wise, but (iii) as x increases still 
further the z's do not follow suit: x = 10^ corresponds to 

z = 500. This suggests that if x is a value assumed by X, 
that Z shares the properties of X in mid-range, i.e. for 
intermediate x-values, but differs from X by having a dispro- 
portionate probability of assuming small values (near zero) , 
or large values (near, but less than, 1/B) . Thus the hazard of 
Z will appear to be a "bath tubbed" version of X, particularly 
if X is exponential. 

We focus attention on the representation (3.6) in what 
follows, mainly for analytical and computational convenience. 

Of course there are many other possibilities, such as 
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L(x) = 1 - e 



-ax 



R(x) 




(3.7) 



these latter may be adjusted to provide sharper-edged bath tubs 
than can (3.6), but iteration of (3.6) may be induced to 
accomplish the same purpose. 

( 2 ) A Decreas ing Fa i lure Rate Mo del 

Define the random variable W in terms of X, X again 
being exponential with parameter X ^ : 

W = XT (X ) (3.8) 



where T(x) is an increasing function of x, L(0) = 1. Then 
the hazard may be made to exhibit a decreasing behavior. 



Example . Suppose 

(3.9) 

(3.10) 



T(x) = 1 + cx 



c > 0 , 0 < x . 



Then 



z = x ( 1 + cx) 



is monotonic, and small x-values lead to comparable z-values 
(especially when c is small) but larger x-values are "amplified" 
by 1 + cx to yield increasingly large z-values. 
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Attention will be focused upon (3.9), although other 



possibilities exist that accomplish the same purpose, namely 
that of lengthening the right tail of the distribution of X 
(simulating outliers, for instance) while leaving the body of 
the distribution virtually unchanged. 

IV. MATHEMATICAL PROPERTIES OF THE "BATH TUB” HAZARD MODEL 

Various analytical properties of the previously described 
models will now be recorded. These provide useful insights into 
the behavior of the random variables Z and the underlying 
(generating) variables X. 

A. Monotonicity; Quantiles 

It is convenient to focus on montonic increasing trans- 
formations, i.e. if 



then in order that the above function be monotonically increasing, 
dz/dx > 0. Observe that logarithmic differentiation of (4.1) 
provides 



z = G(z) = z<p (z) 



(4.1) 




(4.2) 



and thus dz/dx >0 if and only if 
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(4.3) 



1 + 4> ' (x) 
X (p (x) 



> 0 



Alternatively, the condition is, in terms of L(x) and R(x), 



1 , L' (x) R' (x ) 
x L (x) R(x) 



> 0 



(4.4) 



It is easily seen that the important example (3.6) , 

, , , ax 1 

^ x 1 + ax 1 + Bx ' 

yields a monotonic relationship between z and x. The fact 
that this transformation can be easily and explicitly inverted 
(solved for x in terms of z) will be exploited subsequently. 

Of course if z(x) is monotonically increasing then so 
is x(z), the inverse function. The events (Z < z) and 
(X _< x(z)) are equivalent, and so 

P{Z < z} = P{X < x (z) } , (4.5) 



from which it follows that if x^ = x (p) is the p-100% 
of X, i .e . 



P { X < x(p)} = p , 



then 



P { Z < z (p) } = P { Z < z (x (p) ) } = p 



quantile 

(4.6) 

(4.7) 
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and so z (p) , the p*100% quantile of Z is simply obtained from 



z(p) = x(p) <J>(x(p)) = x(p) L (x (p) ) R (x (p) ) (4.8) 

In other words we very easily translate from (points on) the 
inverse distribution o f X to the inverse distribution of Z. 
Explicit representation of the distribution of Z is however, 
not often easily possible. 



B. Hazard and Density Function Relationships 

In order to investigate the relationship between the 
hazards of Z and X, begin by writing 



x (p) 

p = F (x(p)) = 1 - exp [- / h (u) du] 

0 x 



or 

x (p) 

/ h (u)du = - £n(l-p) 
0 x 



Now differentiate with respect to p to find 



(4.9) 



h x ,x(p), SglEl- ji 



(4.10) 



or 



h x ,x(p)) - affey- ' i^p = V x(p,) • its : < 4 - 11 > 



nere h^ and f^ are the hazard and density functions of the 
r.v. X. The relationship (4.11) holds for any distribution, of 
course . 
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Differentiation of (3.5) reveals the connection between 

h and h . From (4.2) 
z x 



*<p> 



i . r (x (p) ) 
x (p) 4 > (x (p) ) 



dx (p) 
dp 



(4.12) 



= [<j>(x(p)) + x (p) <f>'(x(p))] 

From (4.11), applied now to the z-hazard, there results 



h ( Z ( P )) * h - fiTp S T + x( p j ♦’ ' X( P>»1 - (4 - 13) 

Z X 



so 



h (z(p)) = h (x(p)) —7 7 rp — 7 — t rp 7 7 rr (4.14) 

z' v ^' X ^ 4 >(x(p)) + x(p) <f>’(x(p)) 

Multiplication of both sides by 1-p then shows, in view of 
(4.11), that the density functions are similarly related: 

f z (z(p)) = f x (x(p)) ^(p) ) + x'Jp) p ( x ( p ) ) 



Example . 

X is exponential (X) . Then 



h z ( V 4>(x (p) ) + x (p) <f> ' (x (p) ) 



(4.15) 



13 



Now use the specific cj>(x) of (3.6): 



a / \ ax 1 

^ (x) 1 + ax * 1 + 3x 

or, in terms of logarithms, 

Jin <J> (x) = Jin ax - Jln(l + ax) - Jln(l + (3x) , 

so 

<j> 1 (x) _ 1_ _ a _ § _ 1 - agx 

<J> (x) x"l + ax~l + 3x _ x(l + ax) (1 + 6x) ' 



and 



<p (x) + x<J> ’ (x) 



<Hx) 



2 + (a -t B)x 
(1 + ax) (1 + 6x) 



(4.17) 



finally 



h (z (p) ) 



A(1 + gx(p)) 2 (l + gx (p) ) 2 
ax (p) [2 + (a + 6) x (p) ] 



(4.18) 



Although this expression is not quite explicit, qualitative properties 

of h can be deduced from it. 
z 

(a) If p -»• 0 , x(p) = - Jln(l-p) 4 0, and hence 

V z( p>> ~ tAtpT < <4 - 19 > 

or 

lim x (p) h (z (p) ) = -y- (4.20) 

p •*» 0 Z Za 
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Since for p -*■ 0 , 





2 

z (p) ~ ax (p) 


and hence 


x(p) ~ [ z (p) /a] (4.21) 


there results 


h ( Z (p) ) - yjj 

2[z(p)] 1/Z & 


or 


lim /z (p) h (z(p)) = — - — (4.22) 

p -*■ 0 z 2 /a 


This shows that 


h z (z) t oo as z 0 , creating the left-hand end 



of the bath tub of Figure 1. 



(b) If p + 1, 


x(p) + and z (p) + 1/6 so 


or 


2 2 

, , » v , a6 x (p) 

h z (z(p)) ~ X - (a + B) 

1 R 2 

lim - 1 ~ h (z(p)) - A ^ R (4.23) 

p -»• °° [ x ( p ) ] 


For p •> 1 


1 ~ 6z(p, ~a$ + x<p) (4 - 24) 


so 


*<P> l -tfilpT <4 ' 24) 


and thus 
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h (z(p)) ~ A 2_±JL i ^ (4.26) 

Once again it appears that the hazard rises rapidly, this time 
as x(p) + «> and z(p) t 3 the other end of the bath tub 
is thus fashioned. 

(c) If p = 1 - e - '*', then x(p) = A \ Then 

yzu-e- 1 )) - ill ± I 1 + B/Xl 2 (4.27) 

j [2 + (a + 3)/A] 

The bath tub effect is presumably achieved by choosing a large 
and 3 small. Let a 0 and 3 0 independently in (4.27); 

it is clear that the limiting value of the hazard is A. This 
indicates that the hazard is (approximately) A for middling 
values of z . 



C. An Explicit Formula for a Hazard 

The expression (3.6) leads to the relationship 



z (p) 



T T 



2 , 4 

ax (p) 

ax (p) ) (1 + 3x (p) ) 



(4.28) 



and the latter may be explicitly inverted by solving a quadratic 
equation. The result is 
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X (p) 



(g + B) z ( p ) + /(g + B) 2 z 2 (p) + 4z (p) g(l-gz(p)) 

2g (1 - 6z (p) ) 



(4.29) 



Now a direct differentiation of this expression and invocation 
of (4.11) produces the expression 



This form, while explicit, provides no particularly useful insights; 
the bath tub end shapes already noted in (4.22) and (4.26) can be 
deduced directly from (4.30). 



illustrate the behavior of the present hazard representation in 
a more understandable fashion than does the formula itself. 

D. An Explicit Formula for the Failure Time Distribution 

Because x and z are monotonically related through 
(4.28) we have 





+ (g + 6) 




(4.30) 



Some graphical plots of h^ are presented later. They 
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P{Z <z}=P {x<x = ( a+ ^) 2 + /(a+3)^z 2 + 4az(l-3z) 



2a (1 - 3z) 



= 1 - exp { -A 



(a+3) z 



+ *4 



2 2 

a+3 ) z + 4az(l-3z) 



2a (1 - Bz) 



Jl 



(4 .31) 



Again the explicit formula seems unproductive of insights, 



V. MATHEMATICAL PROPERTIES OF THE DECREASING FAILURE RATE MODEL 
A. The Hazard Behavior. 

The expression (4.14) can be applied to deduce the hazard 
function of the representation (3.8), advertised to produce a 
decreasing failure rate. There we specified 

cf>(x) = T(x) = 1 + cx , (5.1) 

and thus, from (4.14), 

h w (w(p)) (1 + cx (p) ) + x ( p) c 1 + 2cx(p) (5.2) 

Qualitative properties follow easily. 

(a) If p ->■ 0 , x(p) 1 0, and 

h (w (p) ) ~ A (5.3) 

w 

Thus the hazard is approximately A for small z. 
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(b) if p 



1, x(p) t <*>, and 



w (p) ~ 



c [x(p) ] 2 



so 

h (w (p) ) ~ — A ■- ■ , (5.5) 

2/cw (p) 



which clearly decreases, as claimed. It may be inferred that the 
distribution of W appears nearly exponential, but possesses 
an extraordinarily long right tail — possibly the result of outliers. 



B. Explicit Formulas for the Hazard and the Distribution Function 
Direct solution of the quadratic equation 



2 

w = x(l + cx) = x + cx 



presents 



x(p) 



/I + 4 cw ( p ) - 1 



which, when differentiated, leads to 



h (w) = h (x) 

w /r + " 4cw x 



/I + 4cw 



(5.6) 



(5.7) 
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The distribution function is 



P{W < w} = P 



X < x 



/I + 4cw - 1 
2c 



(5.8) 



1 - exp < - x 



/I + 4cw - 1 
2c 




This distribution bears a close family resemblance to the Weibull 
distribution 1 - F(w) = exp{-k/w}, especially for large (right 
tail) values of w. 

VI. AN ALTERNATIVE "BATHTUB” HAZARD REPRESENTATION 

The simple parametric model (3.6) leading to a bathtub-shaped 
hazard is by no means the only possibility. We next describe 
another simple approach. It is that of defining a hazard function 
having an appropriate shape, and then deducing the corresponding 
distribution function, and a procedure for sampling from it, 
rather than proceeding in reverse order, as before. 



where g(z) >0 is a decreasing function of z such that 
lim^ ^ ^ g(z) = 0, and k(x) is an increasing function of x, 
such that (preferably) k(0) = 0 and k (<*>) = °». Such a function 
can yield a bathtubbed hazard. 



Let the hazard be of the form 



h(z) = g(z) + X + k(z) 



( 6 . 1 ) 
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Example . 



h (z) = 



z + a 



+ Bz + A 



( 6 . 2 ) 



A, B, a, A all positive. 

Clearly, (6.2) has a generally "bath tub-like" appearance, 

since 



for if 



h* (z) = - 



- - + B < 0, 
a ' 



+ B 



(z + a ) 

then h’ (0) < 0 



(6.3) 



(6.4) 



while for 



z > z. 



■dr 



- a 



(6.5) 



h(z) > 0. 

Detailed behavior is adjustable by choice of the parameters. 
Now the distribution function of time to failure, Z, is obtained 
from (6.2) ; 



P{Z > z} = exp|- 



/ h(x)dx 
0 



= exp ] - / (-■ + ■■ — + Bx + A) dx 



= exp 



H> 



An(l + -) + t z 2 + Az 
a 2 



/ ct i , B 2 . -Az 

- ( stt) exp( ‘ 2 2 > e 



= F 1 (z) F 2 (z) F 3 (z) , 



( 6 . 6 ) 
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where 



A 



F 1 (z) 



a 



a + z 



F 2 (z) = exp (- | z 2 ) 



(6.7) 



F 3 (z) = e 



All of the above are recognized as being the complements of 
distribution functions. In effect, the distribution of Z 
is that of the minimum of three independent random variables: 



having the distribution F^ = 1 - F^ (i = 1,2,3) . This fact 
leads directly to an easy procedure for simulation of Z by 
simply obtaining the smallest from among the realization of 
X^, X 2 , and X^. The advantage of the previous method, based on 
(3.6) for instance, is that only one realization — that of an 
exponential in that specific example — leads to the realization 
of Z. This is not only computationally attractive, but seems 
to facilitate the application of such Monte Carlo variance 
reduction techniques as control and antithetic variables 
cf . Hammersley and Handscomb (1964) . 



P{Z > z} = P{X 1 > z}*P{X 2 > z}*P{X 3 > z) 



( 6 . 8 ) 
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VII. OBTAINING SPECIFIED HAZARD BEHAVIOR BY SIMPLE SAMPLING 



The development of the last section illustrates one manner 
in which hazard behavior may be conveniently represented and 
simulated. We now show how such behavior may alternatively be 
obtained by simple simulation, i.e. from one realization of a 
basic (possibly exponential) random variable. 

Refer to (3.5), in which 



Z = G (X) 

and, if G(*) is monotonically increasing, 



(7.1) 



z (p) = G(x(p) ) , (7.2) 

z (p) and x(p) being the p*100% percentiles of Z and X 
respectively. Then the counterpart to (4.14) that results from 
differentiation of (7.2) is the expression 



h z U(p)> = h x (x(p) 1 (TTxTpTT ° h x (x(p)) (dz/dx) (7 - 3) 



Consequently, if one specifies h z (z) as a suitable function of 

the "time" z, and specifies the distribution of the stochastic 

variable X--and hence its hazard, h — there results a differential 

x 



equation for z (x) = G(x): 



h (z) = h (x) ; 

z dx x 



(7.4) 
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integration then provides the desired transformation, G. In 
other words, we seek z (x) satisfying 



z 



X 



/ h (u) du = / h (v )dv 
0 z 0 x 



(7.5) 



which can sometimes be carried out in a useful closed form. 

Example 7.1 . Refer to the example of Section VI, wherein h 
is given by the expression (6.2) and we assume that X is 
exponential, so h^ is constant. Then in order to determine 
G(x) = z (x) , solve the equation 



Closed-form solution of this expression for z in terms of x 
is of course impossible. One possible approach is purely 
numerical: find an approximate solution, z^ (x) , e.g. the appro- 

priate solution of the quadratic 




or 



A £n ( 1 + — ) + ^ z^ + Az = x 
a 2 



(7.6) 



B 2 . 

2 z + Az = x 



(7.7) 
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and then correct the result by a few Newton-Raphson iterations. 
In other words, put 



z 1 (x) 




(7.8) 



B 



now apply Newton to obtain an improved solution 

A £n(l + z /a) 

z 2 (x) z i a/(z 1 + a) + Bz^ + A (7.9) 

which will be feasible if 0 < z 2 • The process can be iterated 
(the numerator will change after the first iteration) . If one 
wishes to use this model it may actually be desirable to start by 
solving 



for z , in which case the numerator will not be as shown in (7.9) 
convergence may be more rapid. 

Example 7.2 . Change the hazard representation of the previous 
example as follows: let 



§■ z^ + (A + —) z - x 
2 a 



0 



(7.10) 



h 2 (v) 



(v + a) 



— 2 + Bv + A; (A, a, B, A > 0) 



(7.11) 



then 




B 2 

+ 2 z + Az 



(7.12) 
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Now the solution of the equation 



Az 

a ( z + a) 






+ Az = x 



(7.13) 



can be carried out in closed (if cumbersome) form, since (7.13) 
is a cubic. Hence an explicit representation of z (p) in terms 
of an exponential's x(p) can be provided. Perhaps more importantly, 
a direct, simple, simulation of a bathtubbed random variable is at 
hand. 



VIII. CONCLUSION 

The purpose of this report has been to show that rather 
complex distributional behavior can be represented in terms of 
simple standard distributions. We have concentrated here upon 
alterations of the exponential, but it is obvious that other 
distributions can be treated similarly. 

In a later report we will illustrate the use of this technique 
in data analysis, and in the simulation of non-homogeneous 
stochastic point processes. 
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